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Abstract 

Using results from conformal field theory, we compute several universal am- 
plitude ratios for the two-dimensional Ising model at criticality on a symmetric 
torus. These include the correlation-length ratio x* = lim^oo and the 

first four magnetization moment ratios Vi n = (A4 2n ) / (M 2 ) n . As a corollary 
we get the first four renormalized 2n-point coupling constants for the mass- 
less theory on a symmetric torus, G\ n . We confirm these predictions by a 
high-precision Monte Carlo simulation. 
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1 Introduction 



A central concept in the theory of critical phenomena is the idea of universal- 
ity, which states that phase-transition systems can be divided into a relatively small 
number of "universality classes" (determined primarily by the system's spatial dimen- 
sionality and the symmetries of its order parameter) within which certain features of 
critical behavior are universal. In the 1950s and 1960s it came to be understood that 
critical exponents are universal in this sense |IJ. Later, in the 1970s, it was learned 
that certain dimensionless ratios of critical amplitudes are also universal 0. 

The past quarter-century has seen enormous progress in the determination of 
critical exponents for a wide variety of universality classes, including exact analytical 
results for two-dimensional (2D) models || ||, [5], || and increasingly precise numerical 
determinations for three-dimensional models by a variety of techniques (field-theoretic 
renormalization group [0, ||, series extrapolation || [L0|, [11], [L^, [13], [L4], [15], [RJ, Monte 
Carlo |T7|, 0, [DJ |2D], ^ f^, (£| [2§ (§, fZfJ). As a result, attention has turned quite 
naturally to universal amplitude ratios: these include amplitude ratios in infinite 
volume and those in finite-size scaling (FSS). Though much numerical work has been 
done, few exact results are known.[] 

The critical behavior of many 2D models can be studied analytically using con- 
formal field theory (CFT) [|], |5|, ||. Many critical exponents have been determined 

lot 



exactly, along with a few universal amplitude ratios |52], |54|, |56], ^7], |58 



|B"T] , O E3[ The main goal of the present paper is to compute, using CFT, a few 
more universal amplitude ratios for the 2D Ising model and to test these predictions 
by a high-precision Monte Carlo study. The amplitude ratios considered here arise 
in finite-size scaling; they can be computed starting from the correlation functions of 
the critical 2D Ising model on a torus. 

The first class of quantities we study concern the shape of the magnetization 
distribution p(A4) at criticality on a symmetric torus [L x = L y ).f\ We study the 
rescaled shape of this distribution (i.e. normalizing by its width (Ai 2 ) 1 ^ 2 ) as well as 
the dimensionless ratios of its moments, 



2/i 



(M 2n ) 
(M 2 ) n 

We can also define the dimensionless cumulants 

(M 2n ) com 



u 2n = 



(M 2 y 



XV, 



(1.2) 



1 Among the models studied are the 2D Ising model 

~_3§ 




2D Potts models 



B PI 

IIP 



self-avoiding walks |L9|, the 3D Ising model |40|, 41 
3D O(N) spin models |7] ||, g ||, ^6] |g [50| 
^III . This list of references is far from exhaustive 



|28|, |2g |3C|, ^lj, 2D nonlinear cr-models 
the Baxter 8-vertcx model |p9(l, 2D and 3D 
^,P,p]p7]p] P, |2^ 



3D site percolation 



and the 5D Ising model 



2 One of the insights of conformal field theory is that universal finite-size-scaling properties (such 
as the universal amplitude ratios considered here) ought to be studied as analytic functions of the 
modular parameter t of the torus. Nevertheless, we think that the case of a symmetric torus (r = i) 
is of sufficient practical importance in Monte Carlo simulations to warrant special attention. 
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For any symmetric distribution p(M) = p(—M) these satisfy^ 

U 4 = V4-3 (1.3a) 

U 6 = V Q - 151/4 + 30 (1.3b) 

U 8 = V 8 - 28V 6 - 35V? + 42OV4 - 630 (1.3c) 

U w = V w - 45Vb - 2IOV4V6 + 1260V 6 + 3150V; 2 - 18900V 4 + 22680 (1.3d) 



Note that V4 and U4 are closely related to the so-called Binder cumulant [64 



u -1 {M ) -1 ^ ri4^ 

^Binder = 1 " - 1- y - -y . (1.4) 

For (3 < (3 C the ratios V 2n tend in the infinite-volume limit to those characteristic of 
a Gaussian distribution, 

^(Gaussian) = (2n-l)!! (1.5a) 
U 2n ( Gaussian) = (1.5b) 

while for (3 > (3 C they tend to those characteristic of a sum of two delta functions, 

V2n(two deltas) = 1 (1.6a) 

U 2n (two deltas) = - B 2n (1.6b) 

n 

where B 2n = (— l) n_1 (2ra)! £(2n)/2 2n ~ 1 ir 2n is a Bernoulli number. At @ — fi c , however, 
these ratios acquire non-trivial values in-between ( 1.5a ) and (|1.6a|) .H These values, 



which are universal, can in principle be computed by integrating the spin correlators 
for the critical 2D Ising model on a torus, which were determined by Di Francesco et 
al. [^3], [54| using CFT. In practice, however, the formula for V 2n rapidly gets more 
complicated as n grows. Di Francesco et al. [|3|, |54j computed V 4 to roughly three 



decimal places by Monte Carlo integration. Here we shall improve this result by three 
orders of magnitude, and shall also compute V 6 to five decimal places, V 8 to almost 
four decimal places, and Vio to three decimal places: 

V 4 = 1.1679229 ±0.0000047 (1.7) 

V 6 = 1.4556491 ± 0.0000072 (1.8) 

V 8 = 1.89252 ±0.00018 (1.9) 

Vio = 2.53956 ±0.00034 (1.10) 



3 These relations can be computed from the generating functions 

C/gn 

(2n)V 



OC jj I oo 




n=l 

with Vo = V% = 1 and V2n+i = 0. 

4 The Schwarz inequality implies that V% n > 1 for any model, and the Gaussian inequality |]65| , |6r| 
implies that Vi n < (2n — 1)!! for ferromagnetic Ising models. In particular, we have —2 < C/4 < 
0. Moreover, Newman [B7| and Shlosman [pSfl have proven, for ferromagnetic Ising models, that 
{-l) n - 1 U 2 „ > for all n; and Newman |67) has proven some additional inequalities on the Ui n - 
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Finally, we shall measure the ratios V2 n for 2 < n < 10 by Monte Carlo simulation, 
with an accuracy that gradually deteriorates as n grows. For n = 2, 3, 4, 5 our Monte 
Carlo estimates agree well with the theoretical predictions (but are of course less 
precise). 

Another interesting quantity is the second-moment correlation length £. It has a 
sensible definition in finite volume (see Section ||), and its expected FSS behavior is 



where the leading coefficient 



x 



x* + AL~^ + 



hm ijL 

L— too 



;i.ii) 



;i.i2) 



is universal. J 5 } f] Here we shall compute x* for the 2D Ising model at criticality by 

II; we find 



numerically integrating the known spin correlators [53 



x 



0.9050488292 ± 0.0000000004 



(1.13) 



Our Monte Carlo data confirm this prediction. To our knowledge, this is the first 
exact determination of x* for any universality class. Monte Carlo estimates of x* are 
available for many other 2D models — including the 3-state Potts model |75 |, the 



4-state Potts model [76], f77[ , the 3-state square-lattice Potts antiferromagnet [78], [79 
the XY model j8C§ , and several points on the self-dual curve of the symmetric Ashkin- 
and it would be very interesting to compute x* analytically for 
Numerical estimates of x* are also available for some three- 



Teller model [76 
some of these models, 
dimensional spin models 



18, 41, 47, 24 



Consider, finally, the dimensionless renormalized four-point coupling constant 

U 4 



94 



U4 



(£/£) c 



1.14) 



where u<± is the connected four-point function at zero momentum. In the FSS limit 
L — > oo, (3 — > f3 c with £/L fixed, g^ is a nontrivial function of the FSS variable 



g A = F 9A {i/L) 



;i.l5) 



Therefore, the function g^.L) fails to be jointly continuous at ((3, L) = (j3 c , oo); 
many limiting values are possible depending on the mode of approach, and the massive 
and massless scaling limits 



9l = il 1 ? r lim 9i(P, L) 
G l = r lim 1 ™94(P,L) 



lim g<Sfic, L) 



;i.l6) 
;i.l7) 



5 The quantity x* also plays an important role in a recently developed method for extrapolating 
finite- volume Monte Carlo data to infinite volume g [70| ^ ||, f72|. 

6 An analogous situation holds in a cylindrical (L x oo) geometry for the exponential correlation 
length in the longitudinal direction, £, e xp(L) [which can be defined in terms of the logarithm of the 
ratio of the two largest eigenvalues of the transfer matrix] . Privman and Fisher |73|| showed that 

lim £,exp{L)/L at criticality is universal, and Cardy |74| showed that for 2D conformal- invariant 

L^oo — 

systems it is equal to 1/(^77). 
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correspond to the two extreme cases gl = F g4 (0), Gl = F 9i (x*). As a corollary of 
our computation of V4 and x*, we obtain the value of g^ at criticality on a symmetric 
torus: 

d = --pi = 2.2366587 ±0.0000057. (1.18) 

More generally, consider the dimensionless renormalized 2n-point coupling con- 
stant 

= ^ri^ , (1-19) 

where r 2n is the amputated one-particle-irreducible 2n-point function at zero momen- 
tum.] 7 ] We can predict the next three renormalized coupling constants at criticality 
on a symmetric torus: 

u _ in/7 2 

G* 6 = -— — — = 29.25457 ±0.00015 (1.20a) 



x 



Gg = - V > ~ ^ ± 280U * = 942.6095 ±0.0072 (1.20b) 



G 



U w - 120U 4 U 8 - 126C/| + 4620£/ 4 2 £/ 6 - 15400t/ 4 4 



10 " x * 8 



56110.24 ±0.56 (1.20c) 



In addition, we shall provide Monte Carlo estimates of G\ 2 through [cf. ( [4.20D 
below] . 

This paper is organized as follows: In Section ^| we review the relevant exact 
results available for the 2D Ising model at criticality on a torus, and we compute (by 
numerical integration) the CFT prediction for the quantities x*, V±, V 6 , V 8 and V w . In 
Section |3| we explain the Monte Carlo algorithm we have used to simulate this model. 
In Section |] we analyze our numerical results for the static observables and compare 
them against the available exact results. Finally, in Section [5] we present our final 
conclusions and discuss prospects for future work. In Appendix |A] we explain how we 
carried out the numerical integrations involved in computing x*, and in Appendix [FJ 
we summarize the definitions and principal properties of the Jacobi theta functions. 



2 Theoretical Results 



The universal amplitudes we consider in this paper (x* and V^) can be written in 
terms of integrals of the 2n-point spin correlation functions of the the critical 2D Ising 



7 The are defined by the generating-function relation 

f 



E 



2n ,2n 



n=l y ' 

Recall that u 2 — X an d that u 2n = L~ d (M 2n 



sup 

J 



U 2n j2n 



1 (2n) 
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continuum field theory on a torus. These correlators were obtained by Di Francesco 
et al. |)3| [54|] using an approach based on conformal field theory. The result is 8 



YH=i Z u (a Zl ■ ■ ■ (J Z2n )v 



where 



2\r]\ 



and 



'- I ' °~Z2n ) V 



1 y 

On+21^12 Z-^ 
z I '/I £ ,=±i 



£^=o 



^(o)^ 2 

2 n + 2 |r/| 2 



E 

e j =±i 



0„ 



n 



6»i (zj - 2,- 



^i(O) 



Eitj/2 



i<j 



\HiHjl2 



(2-1) 



(2.2) 



(2.3a) 



(2.3b) 



Here we have used the complex-number notation z = xi + 2x2; #i(0) ~ 2.8486946040 
is the derivative of t) with respect to z evaluated at z = and r = i; and 
77 ~ 0.7682254223 is the usual Dedekind function r/(r) evaluated at r = i. Please 
note that ^(0) = 2mf [cf. ([BT3|) 1. Note also that the contribution of {ej} to (2.3) 
is equal to that of {— Sj}, so in the numerical evaluation of this expression we need 
only take half the terms. The expression ( |2.1| ) gives the FSS limit for the Ising-model 
correlation functions at criticality: here Zi denotes the position in lattice units divided 
by the lattice linear size L. 

Remark. Although the sector u — 1 does not contribute to the partition func- 
tion [since Z\ ~ #i(0) = 0], it does contribute to the correlation functions [since 

z\ '" <J z 2n )i 7^ 0]. So this sector cannot simply be discarded. See ref. 
details. 



53 for 



The two correlators that are needed in the evaluation of the Binder cumulant are 



Z v {o Zl cy Z2 ) u 





)) 


1/4 










2 


V 




0i (z% - z 2 ) 


1/4 



(2.4) 



Z\ &Z2 &Z's &Z4 1 1 



2V2I77I 



1/2 



Z\ + Z 2 — Z3 — Z4 



8 There is a misprint in the normalization of the 4-spin correlator in equation (9) of |p4| , and in the 
normalization of the 2n-spin correlator in equation (6.6) of |5q| . We have rederived both correlators 
using the chiral bosonization prescription presented in |p3[ . With the correct normalization, shown 
in (2A)-(2.3) below, we are able to reproduce the numerical value of V4 reported in [^ij, as well as 
the numerical estimates of V4, Ve, Vs and Vto obtained in our simulation. 
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0i(zi - z 2 )9 1 (z 3 - z 4 ) 



1/2 



0i (21 - 23)01(21 - 2 4 )0i(z 2 - 23)01(22 - Z4) 

x 1/2 

+ (2 <-> 3) + (2 <-> 4) (2.5) 



We also need the 6-point correlator to compute V 6 . Its exact expression can be 
deduced easily from the general equation (2.3): 



Z„(a Zl a Z2 a Z3 a Z4 a Zs a Z6 ) l 



K(Q)I 3/4 
4(77 1 



0, 



2l + Z2 + Z3 — Zi — 25 — Zq 



x (z ls za, z 3 , z 4 , z 5 , ze) + (2 <-> 4) + (2 <-> 5) + (2 <-> 6) 

+ (3 <-> 4) + (3 «-> 5) + (3 «-> 6) + (2 «-> 4; 3 <-> 5) 

1 1/2 

+ (2 «-> 4; 3 «-> 6) + (2 «-> 5; 3 «-> 6) [ 



where the function \1/ is defined as 



, Z 2 , Z3, Z4, Z5, Zq) 



'01(212)01(213)01(223)01(245)01(246)01(256) ' 



1/2 



2,3;i=4,5,6 @l( Z ij) 



(2.6) 



(2.7) 



and we have used the shorthand notation = Zi — Zj. 

From these equations we can obtain the values of x* = hin^oo £/L and V% n by 
numerical integration. In particular, the correlation length on a periodic lattice of 
size L is defined to be 



4-1 



1/2 



(2.8) 



2sin(7r/L) \F 

where \ is the susceptibility (i.e., the Fourier-transformed two-point correlation func- 
tion at zero momentum) and F is the corresponding quantity at the smallest nonzero 
momentum (2vr/L,0) [see (3.7)/(3.9) and (gH/(B) below for details]. This is 
a finite-lattice generalization of the second-moment correlation length. Then, the 
universal amplitude x* is given by 

vV2 



x 



where 



X 
F 



X 

2tt \F ' 

d 2 z (a a z ) 

d 2 z (a a z ) cos(27rxi) 



(2.9) 

(2.10) 
(2.11) 



and J d 2 z = f J dx\ dx 2 - The details of this computation are given in Appendix 
We obtain 



d 2 z (a a z ) 



1.55243295465 ± 0.00000000004 



(2.12) 



J d 2 z (a a z ) cos(2ttxi) = 0.04656744682 ± 0.00000000004 (2.13) 



As a result, we obtain x* with 10 digits of precision: 



x * = 0.9050488292 ± 0.0000000004 . (2.14) 

We repeated the computation requiring 11 digits of precision in the integrals, and the 
result was the same. 

The universal moment ratio V4 is given by 

v ^ _ / d 2 z 2 d 2 z 3 d 2 z 4 {a a Z2 a Z3 a Z4 ) ^ ^ 

[J d 2 z (a a z )} 2 



Di Francesco et al. [53, 54] performed the integrals in numerator and denominator 
by Monte Carlo and obtained 

V A = 1.168 ±0.005 . (2.16) 



We have improved this value, as follows: For the denominator of ( 2.15 ), we use the 
very precise estimate Q2.12D coming from deterministic numerical integration. For the 
numerator, we performed a Monte Carlo integration using 10 9 measurements. Our 
result is 

V 4 = 1.1679229 ±0.0000047, (2.17) 

which is compatible with (|2.16|) but three orders of magnitude more precise. This 
value also agrees closely with the estimate of Kamieniarz and Blote based on 
extrapolation of the exact results (computed by transfer-matrix methods) for L < 17: 

V A = 1.1679296 ± 0.0000014 , (2.18) 

where the error bar is of course somewhat subjective .0 

More generally, the universal moment ratio V 2n is given by 

T/ Jd 2 z 2 ■■■ d 2 z 2n (a a Z2 ---(T Z2n ) 

V2n - rr „ , Tm • 

[J d l Z {<7 (7 Z )\ 

We have been able to compute the (exact except for the numerical integration) values 
of the ratios V 6 , V 8 and V w . We performed the integrals in the numerator by Monte 
Carlo, using 10 9 measurements for Vq, 4 x 10 6 measurements for V§ and 2.5 x 10 6 
measurements for Vio- We obtain 

V 6 = 1.4556491 ±0.0000072 (2.20) 
V 8 = 1.89252 ±0.00018 (2.21) 
Vio = 2.53956 ± 0.00034 (2.22) 

In general, the formula for the 2ra-point function contains (2n)!/[2(n!) 2 ] terms [this 
takes into account the {e^} <-» {— ej} symmetry], and this grows asymptotically like 

9 Unfortunately, Kamieniarz and Blote [^l) reported only meager details of the fits that led to 
this extraordinarily precise estimate. That is a shame, as information on the presence or absence of 
particular correction-to-scaling terms could be of considerable theoretical interest. 
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4". Thus, in computing V4 (resp. Vq, Vg, Via) we had to include 3 (resp. 10, 35, 126) 
terms, and the computation of V\ 2 would require handling 462 terms. Moreover, the 
numerator has to be integrated over a (An — 2)-dimensional torus. These facts make 
the high-precision numerical integration of V 2n extremely time-consuming as soon as 
n becomes moderately large. 

Let us consider, finally, the dimensionless renormalized four-point coupling con- 
stant #4 defined by 

where U4 is the connected four-point function at zero momentum, and more generally 
the dimensionless renormalized 2n-point coupling constant g 2n defined by 



92r. 



■in 



(2.24) 



where T 2n is the amputated one-particle-irreducible 2n-point function at zero mo- 
mentum. In the FSS limit L — > 00, (3 — > (3 C with £/L fixed, g 2n becomes a nontrivial 
function of the FSS variable £/L, 



92n = F g2n (£/L) . 



(2.25) 



(There is some evidence that F gi is a decreasing function of In particular, the 

massive and massless scaling limits 



9*2n = fe? r Um 92ll(P) L) 
G*2n = r lim ^92n(P,L) 



lim g 2n (Pc,L) 

L— >oo 



correspond to the two extreme cases g 2n = F g2n (0), G 2n 
rently available estimates for the 2D Ising model are 



F 92n (x* 



(2.26) 
(2.27) 

The best cur- 



9i 



' 14.694 ± 0.002 by high-temperature expansion [84, 28 
14.66 ± 0.42 by e-expansion |28] 
15.50 ± 0.84 by (^-expansion [55 



(2.28) 



14.66 ± 0.06 by expansion around d = [55, |S7| 
I 14.7 ±0.2 



by Monte Carlo p 



Lt I 



2.239 ± 0.007 by Monte Carlo j$| 



(2.29) 



10 Baker and Kawashima [Q conjecture that g4,((3,L) [resp. £(/3, L)\ is a decreasing [resp. in- 
creasing] function of (3 for each fixed L < 00; these two facts, if true, would immediately imply 
that -F 94 (£/L) is a decreasing function of its argument £/L. Numerical data for the 2D Ell Q and 
three-dimensional [T4l], Q Ising models clearly support the Baker-Kawashima conjecture. 
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96 



' 794.1 ± 0.6 by high-temperature expansion [27 , p8[ 
797 ±9 by e-expansion [^8], ^] 



792 ± 40 
691 ± 29 
I 850 ± 25 



by (yf-expansion with constrained g\ |p0 



(2.30) 



by expansion around d = |86|, |87| 
by Monte Carlo 



29.34 ± 0.20 by Monte Carlo |5T| 



(82.5 ± 0.6) x 10 3 by high-temperature expansion |28| , |29 
(83.8 ± 3.2) x 10 3 by e-expansion |28], [291 
(89 ± 5) x 10 3 by Monte Carlo P[| 



#10 



947 ± 10 by Monte Carlo [31 



(12.8 ± 0.7) x 10 6 by high-temperature expansion [28. 29| 
(8.0 ± 1.4) x 10 6 by e-expansion ER |29fl 



(2.31) 

(2.32) 

(2.33) 
(2.34) 



Our own Monte Carlo data, reported in Section [4.3| , combined with the theoretical 
value ( pTBD for x*, improve (|f29|)/(p3lD/(p33|) to G* = 2.23685 ± 0.00016, G* 6 = 
29.2602 ± 0.0047 and G* 8 = 942.91 ± 0.25, respectively, and also give values for G* w 
through G* 20 [see ( ^20|) and the footnote following it]. From ( gig ) and (pTT7|) (gT22f) 
we obtain the theoretical predictions 



i 

6 



= 2.2366587 ±0.0000057 



I7 fi 



10t7? 



29.25457 ±0.00015 



U 8 - 56U 4 U 6 + 280L/f 



.r 



*6 



942.6095 ± 0.0072 



G* 



U 10 - 120U 4 U 8 - 126C/| + 4620C/|C/ 6 - 15400f/ 4 4 



10 



*8 



56110.24 ±0.56 



(2.35a) 
(2.35b) 
(2.35c) 

(2.35d) 



The error bars in ( |2.35| ) are obtained by carefully propagating the statistical errors 
from the Vm [in which the errors are independent except for an extremely small effect 
arising from the value of the denominator ( |2.12j )] to the U2n (in which the errors are 
correlated) and thence to the G* 2n . 
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3 Numerical Simulations 



We consider the two-dimensional nearest-neighbor Ising model on a L x L square 
lattice with periodic boundary conditions, given by the Hamiltonian 

Rising = (3.1a) 

1 (H) 

= -f3 E <W f + const . (3.1b) 

(ij) 

Note that we use throughout this paper a non-standard normalization of (3, which is 
motivated by considering the Ising model as a special case of the g-state Potts model; 
it differs by a factor of 2 from the usual Ising normalization. In our normalization, 
the critical point is at 

(3 C = logtl + v^) « 0.881373587. (3.2) 



3.1 Observables to be measured 

We have performed simulations of this system using the Swendsen- Wang algorithm 
|88| , 89| , p0|| . In particular, we have measured the following basic observables: 

• the energy density (i.e., the number of unsatisfied bonds) 

8 = £(1-<WJ (3.3) 

• the bond occupation 

M ee E n xy (3.4) 

(xy) 



the nearest-neighbor connectivity (which is an energy-like observable ||75|| ) 

where j xy equals 1 if both ends of the bond {xy) belong to the same cluster, and 
otherwise. More generally, the connectivity 7^ can be defined for an arbitrary 
pair i, j of sites: 

/ r \ J 1 if i is connected to j , , 

1tj\\ n ) 1 if i is not connected to j 

the squared magnetization 

= (5>«1 (3.7a) 



X 



q ^{^, V V 2 



1 a =l \ x 



(3.7b) 



11 



where <r x = E R 9 1 is the Potts spin in the hypertetrahedral representa- 
tion^] and V = L 2 is the number of lattice sites 



powers of the squared magnetization 

M 2n = (M 2 ) n 



(3.8) 



the square of the Fourier transform of the spin variable at the smallest allowed 
non-zero momentum 



E CT z e 

X 

1 q 



2mxi/L 



E°"* 



^2irix2/L 



q 



1 2 ^ 



D 2mxi/L 



+ 



^2irix2/L 



(3.9a) 
(3.9b) 



where (xi, 22) are the Cartesian coordinates of point x. Note that T is normal- 
ized to be comparable to its zero- momentum analogue M 2 . 



the mean-square and mean-fourth-power size of the clusters 

S2 = E#(^) 2 

c 

= E#(£) 4 



(3.10) 
(3.11) 



where the sum is over all the clusters C of activated bonds and #(C) is the 
number of sites in the cluster C. 

From these observables we compute the following expectation values: 

• the energy density E per spin 



E 



the specific heat 



Ch = ^var(£) = I [(S 2 ) - (S) 2 



the magnetic susceptibility 



X 



the higher magnetization cumulants 



U2n 



±(M 2 ) 



V (M 2n ) 



(3.12) 

(3.13) 
(3.14) 

(3.15) 



11 Let {e(")}^ =1 be unit vectors in R 9 " 1 satisfying • = (qd afj - l)/{q - 1), and let 
cr x = e^" x \ For q = 2 this means cr x = cos(7rcr x ) = ±1. 
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• the magnetization moment ratios 



•in 



(M 2n ) 
(M 2 ) n 



the correlation function at momentum (2-k/L, 0) 

1 



V 



the second-moment correlation length 

{ 



1/2 



(3.16) 



(3.17) 



X _ 
2sin(7r/L) \F 

the variant second-moment correlation length 

? 2ir \F ) • 
which differs from £ only by correction-to-scaling terms of order L~ 



(3.18) 



(3.19) 



Remarks. 1. Using the Fortuin-Kasteleyn identities |)T], |92], |93], |9C§, it is not difficult 
to show that 



(AO 
(M 2 ) 



p(B-{8)) 



q 



q-l 
(ft) 

g + 1 



9 



(s-(O) 



(5 2 2 ) 



q 



(3.20) 
(3.21) 

(3.22) 
(3.23) 



where p = 1 — and 5 = 2U is the number of bonds in the lattice. As a check on 
the correctness of our simulations, we have tested these identities to high precision, 
in the following way: Instead of comparing directly the left and right sides of each 
equation, which are strongly positively correlated in the Monte Carlo simulation, a 
more sensitive test is to define new observables corresponding to the differences (i.e., 
M — p(B — S) and so forth). Each such observable should have mean zero, and the 
error bars on the sample mean can be estimated using the standard error analysis 
outlined below. The comparison to zero yields the following x 2 values: 



For (3.20) 
For (fm]) 
For ( gjp 

For rra 



X 2 = 10.23 (14 DF, level = 75%) 
X 2 = 17.00 (14 DF, level = 26%) 
X 2 = 9.93 (14 DF, level = 77%) 
X 2 = 9.05 (14 DF, level = 83%) 



(3.24) 
(3.25) 
(3.26) 
(3.27) 
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Here DF means the number of degrees of freedom, and "level" means the confidence 
level of the fit (defined at the beginning of Section f| below). The agreement is 
excellent. 

2. We also compared our data for V4 for L = 4, 6, 8, 12, 16 with the exact values 
computed by Kamieniarz and Blote [KI]] using transfer-matrix methods. We get x 2 = 



5.14 (5 DF, level = 40%), indicating good agreement. 

For each observable O discussed above we have measured its autocorrelation func- 
tions in the Swendsen-Wang dynamics, 

Cbo(f) = (0,0, +1 > - (Of (3.28) 

*»(*) = £4t < 3 ' 29 ' 

Oe>e>(Uj 

where the expectations are taken in equilibrium. From these functions we have esti- 
mated the corresponding integrated autocorrelation time 

2 00 

Tint.o = » J2 Poo{t) (3.30a) 

^ t=-oo 
1 00 

= 7; + T^Pooit) (3.30b) 
z t=i 

by the methods of Ref. [[54], Appendix C], using a self-consistent truncation window 
of width 6Tj n t,o- This autocorrelation time is needed to compute the correct error bar 
on the sample mean O. 

Remarks. 1. The error bar of the second-moment correlation length is computed 
by considering the random variable 

M 2 T 

a = — - — , (3.31) 

Pm 2 Pr 

which automatically has zero mean. Then, 

where £ denotes our Monte Carlo estimate of £. In practice, the values of hm 2 an d fir 
are replaced by their corresponding sample means (which should be computed first). 
2. The error bar on the ratio Vm is computed in a similar fashion: 

var(Vk) 1/2 = ^var(^) 1 / 2 , (3.33) 

where V 2n denotes our Monte Carlo estimate of V-m, and the observable 0'{ n is defined 
as 

M 2n M 2 

Q'i n = — n— + n - 1 (3.34) 

pM 2n PM 2 
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and has mean zero. Again, the mean values fiM 2n an d Pm 2 are replaced in practice 
by their sample means. 

3. As a further check on the correctness of our simulations, we have computed 
both sides of the identity 



1 



(l-p)(2-£Q 
P C H + (l-p)(2-E) 



(3.35) 



proven in |95], equation 7] (see also f75fl).p| This is a highly nontrivial test, as it relates 
static quantities (energy and specific heat) to a dynamic quantity (autocorrelation 
function of the bond occupation at time lag 1). We have also checked with great 
accuracy the identities WE 



C ££ (t) 
Pee(t) 



-^CW(* + 1) 

PMM(t + 1) 



Paw(1) 



q 



P£8(t + 

Pss{l) 



C ££ {t + 1) 



(3.36) 
(3.37) 

(3.38) 
(3.39) 



3.2 Summary of the simulations 

We have run our Monte Carlo program on lattices with L ranging from 4 to 512 
(see Table H). In all cases the initial configuration was random, and for L < 64 (resp. 
L > 96) we discarded the first 5 x 10 4 (resp. 10 5 ) iterations to allow the system to 
reach equilibrium; this discard interval is in all cases greater than lO 4 ri nt> £-.0 The 
total number of iterations ranges from 2.15 x 10 6 (L = 4) to 8.2 x 10 6 (L = 512), and 
is selected to be approximately 10 6 Ti n t,£- These statistics allow us to obtain a high 
accuracy in our estimates of the static and dynamic quantities (error < 0.17% and ^ 
0.51%, respectively). The static data are displayed in Table [1] (%, F, £) and Table |2] 
(the ratios V2 n ). The dynamic data will be reported elsewhere. 

The CPU time required by our program is approximately 6.3 L 2 fis per iteration 
on a Linux Pentium machine running at 166 MHz. The total CPU time used in the 
project was approximately 7.5 months on this machine. 

12 Please note that |95| used a definition of energy that is slightly different from the one used here: 
£(Ref. U) - 0-/V)(T, (xy) 8,.,* v ) =2 — E. 

13 Such a discard interval might seem to be much larger than necessary: 10 2 Tj nt would usually be 
more than enough. However, there is always the danger that the longest autocorrelation time in the 
system may be much larger than the longest autocorrelation time that one has measured, because 
one has failed to measure an observable having sufficiently strong overlap with the slowest mode. As 
an undoubtedly overly conservative precaution against the possible (but unlikely) existence of such 
a (vastly) slower mode, we decided to discard up to 2% of the entire run. This discard amounts to 
reducing the accuracy on our final estimates by a mere 1%. 
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We have improved the precision of our analysis of the correlation length £ by 
supplementing our own Monte Carlo data with comparable data from Ballesteros et 
al. p6[ . They performed single-cluster |57J simulations of the 2D site-diluted Ising 
model at various concentrations p. Their data for p = 1 (i.e., the usual Ising model) 
correspond to anywhere from 4x 10 5 to 7x 10 5 statistically independent measurements 
at each lattice size from L = 12 to L = 512 (see Table [$J). The statistical independence 
of two consecutive measurements was achieved by allowing 100 single-cluster moves 
between them. Their error bars are slightly larger than ours. As a matter of fact, 
their error bars and our error bars cr(£) satisfy approximately the relation 



*(0 2r int , ,iV' 

W) = V^v— ' (3 ' 40) 

where N (resp. N') is the number of measurements of our (resp. their) work, and O' 
is the observable ( |3.31|) we used to compute the correct correlation-length error bar. 



This supports the belief that their measurements are indeed essentially independent 
and that their error bars are correctly computed. Comparison of their raw data to 
ours at the eleven overlapping L values yields x 2 — 10.28 (11 DF, level = 51%). The 
two data sets are therefore compatible. The corresponding merged data are shown in 
Table |. 



4 Data Analysis 

For each quantity O, we carry out a variety of fits using the standard weighted 
least-squares method. As a precaution against corrections to scaling, we impose a 
lower cutoff L > L min on the data points admitted in the fit, and we study systemat- 
ically the effects of varying L min on the estimated parameters and on the x 2 value. In 
general, our preferred fit corresponds to the smallest L m i n for which the goodness of 
fit is reasonable (e.g., the confidence leveQ is ^ 10-20%) and for which subsequent 
increases in L min do not cause the x 2 t° drop vastly more than one unit per degree 
of freedom (DF). 



4.1 Corrections to scaling 

In the data analysis we should take into account the effect of corrections to scaling 
in order to get reliable estimates of the physical quantities. In particular, the value 
at criticality of any observable O(L) is typically given for large L by 

0{L) = AL Po (l + A'L~ A + ■■■) (4.1) 

14 "Confidence level" is the probability that x 2 would exceed the observed value, assuming that 
the underlying statistical model is correct. An unusually low confidence level (e.g., less than 5%) 
thus suggests that the underlying statistical model is incorrect — the most likely cause of which 
would be corrections to scaling. 
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where po is the critical exponent associated to the observable O, A is the leading 
correction-to-scaling exponent, and the dots indicate higher-order corrections. 

In finite-size-scaling (FSS) theory |98[ for systems with periodic boundary condi- 
tions, three simplifying assumptions have frequently been made: 

(a) The regular part of the free energy, / reg , is independent of L |58| (except possibly 
for terms that are exponentially small in L). 

(b) The relations connecting the nonlinear scaling fields gt and Qh to the conven- 
tional thermodynamic parameters t = j3 c — (3 and h are independent of L |99| 



(c) The scaling field gi associated to the lattice size equals L 1 exactly, with no 
corrections L~ 2 , L -3 , . 
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Moreover, in the nearest-neighbor spin- 1/2 2D Ising model, it has further been as- 
sumed that 



(d) There are no irrelevant operators ||100| , |101| . 



This latter assumption has been confirmed numerically (in the infinite- volume theory) 
through order t 3 , at least as regards the bulk behavior of the susceptibility ||101|| . How- 
both numerical 



ever. 



102, 103 and theoretical 



10 4|| evidence has recently emerged 
suggesting that irrelevant operators do contribute to the susceptibility at order t A . 

The absence of irrelevant operators implies that the corrections to scaling in this 
model are due to the smooth but in general nonlinear connection between the conven- 
tional thermodynamic parameters t and h and the renormalization-group nonlinear 
scaling fields | |105| , |106| , |100| , [LOlj] . Starting from the FSS Ansatz for the Ising- model 
free energy and using the above assumptions, it is possible to obtain a FSS expres- 
sion for the usual observables at criticality as functions of the lattice size L [107|. In 
particular, the leading correction term in the expansion of the susceptibility is the 
L-independent term coming from the regular part of the free energy. This implies 
that for this observable 

A = I ■ (4.2) 



The same result is plausible for the observable F defined in (|3.17|) ; thus, we expect 
A = 7/4 for the second- moment correlation lengths £ and £' and the corresponding 
amplitude x*. The expansion for the magnetization cumulant U2 n gives an exponent 
A = I+7/1/ = 11/4 (perhaps with a multiplicative logarithmic correction). Thus, we 
expect that the ratios Vi n also have a correction-to-scaling exponent given by (|4.2j ) 



[due to the power of the susceptibility appearing in its definition (|1.1|) 1. For a more 
detailed theoretical and numerical analysis of the corrections to scaling in this model, 



sec 



1071 
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4.2 Second-moment correlation length 

The second-moment correlation length £ and its variant £' [cf. ( [3.18 ) / ( 3.19 )1 are 
expected to behave as 



x* + AL~ A + 



(4.3) 



with p — 1. We can estimate p by ignoring correction-to-scaling terms and performing 
a simple power-law fit. We get 

For £: p = 0.99974 ± 0.00036 (L min = 32, x 2 = 1-48, 6 DF, level = 96%) 

(4.4) 

For p = 1.00018 ± 0.00036 (L min = 32, x 2 = 1-37, 6 DF, level = 97%) 

(4.5) 

The agreement with the theoretical prediction is excellent. 

The value of the constant x* can be estimated most simply by fitting the ratio 
£/L or to a constant, ignoring corrections to scaling. We get 

For f : = 0.90577 ± 0.00028 (L min = 32, x 2 = 2.02, 7 DF, level = 96%) 

(4.6) 

For f: = 0.90557 ± 0.00026 (L min = 16, x 2 = 2.73, 9 DF, level = 97%) 

(4.7) 

The estimate based on £ lies 2.6 standard deviations away from the value x* ~ 0.90505 
predicted by CFT [cf. ( |2.14|) 1. The estimate based on £' is slightly better: it lies two 
standard deviations away from the theoretical prediction, and works also for smaller 
L min . Indeed, the corrections to scaling in £'/L are negligible (compared to our 
statistical error) already for L > 16 (see the last column of Table This fact 

makes it almost hopeless to study corrections to FSS in 

If we fit £ to ( |4.3|) , keeping the first correction-to-scaling term and trying to 
estimate simultaneously the three parameters x*, A and A, a good fit is obtained for 

Lmin 8. 

x* = 0.90546 ± 0.00033 (4.8a) 

A = 0.75 ±0.30 (4.8b) 

A = 1.76 ±0.18 (4.8c) 

with x 2 — 2.97 (9 DF, level = 97%). The value of x* is again two standard devi- 
ations away from the theoretical prediction ( 2.14j) . The estimate of A is very close 
to 7/4, and is only 1.4 standard deviations away from 2; but perhaps this estimate 
ought not be taken too seriously, as the correction-to-scaling amplitude is only 2.5 

15 Table || is based on merging our data with that of Ballesteros et al. Q; but virtually identical 
results are obtained using our data alone. 
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standard deviations away from zero (a deviation that is, moreover, comparable to the 
discrepancy in x*). The analogous fit for £' is even more hopeless (the amplitude A 
is compatible with zero within 0.7 standard deviations), so we omit the details. This 
correction-to-scaling exponent A « 2 can be understood as arising simply from the 
ratio = \{L/ix) sin^/L)]- 1 = 1 + (tt 2 /6)L~ 2 + . . . . Indeed, if we fit the data to 
the Ansatz £/L = x* + AL~ 2 we get for L min = 8: 

x* = 0.90569 ±0.00027 (4.9a) 
A = 1.269 ±0.064 (4.9b) 

with x 2 = 4.56 (10 DF, level = 92%). Then, A/x* « 1.40, which is not far from 
tt 2 /6 « 1.64. 

We can improve the precision of our numerical estimates by using the merged data 
of Table ^ (our data plus that of Ballesteros et al. f96|). The simple power- law fit 
yields 

For f : p = 1.00047 ± 0.00033 (L min = 48, \ 2 = 0.86, 5 DF, level = 97%) 

(4.10) 

Forf: p= 1.00074 ±0.00033 (L min = 48, x 2 = 0-82, 5 DF, level = 98%) 

(4.11) 

The fits to a constant give 

For f : = 0.90565 ± 0.00023 (L min = 48, \ 2 = 2.92, 6 DF, level = 82%) 

(4.12) 

For = 0.90555 ± 0.00020 (L min = 16, x 2 = 6.99, 9 DF, level = 64%) 

(4.13) 

The three-parameter fit £/L = x* ± AL~ A is good for L min = 8: 

x * = 0.90552 ±0.00026 (4.14a) 
A = 0.86 ±0.29 (4.14b) 
A = 1.82 ±0.15 (4.14c) 

with x 2 — 7.57 (9 DF, level = 58%). The analogous fit with £' yields a correction-to- 
scaling amplitude compatible with zero within errors. 

In conclusion, one can extract accurate estimates of the critical exponent p and 
the amplitude x* using our Monte Carlo data; the results agree with the theoret- 
ical prediction ( p. 14 ) within two standard deviations. However, it is very difficult 



to estimate from our numerical data the correction-to-scaling exponent (or the corre- 
sponding amplitude). Indeed, for £' the corrections to scaling are negligible (compared 
to our statistical error) for L > 16. 

4.3 Magnetization moment ratios 

If we study the magnetization distribution p{M) as L — > oo at fixed j3, we expect 
three distinct behaviors depending on the value of (3: 
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(a) At (3 < (3 Cl we are in the high-temperature regime, where correlations decay 
exponentially. A variant of the central limit theorem ||108|| guarantees that the 



finite-L distributions will converge, after rescaling by the factor y/Vx, to a 
Gaussian distribution of mean zero and unit variance.^] 

(b) At (3 > P c we are in the low-temperature regime, and the finite-L distributions 
should converge, after rescaling by the factor VMq (where Mq is the sponta- 
neous magnetization), to the sum of two delta functions. There are Gaussian 
fluctuations around these two delta functions, but their width is much smaller, 
namely y/Vxo, where xo is the susceptibility in a pure phase. 

(c) At (5 — (5 C [or more generally, at fixed value of the FSS variable L l l u (f3 — (3 C )], 
the finite-L distributions will converge, after rescaling by the factor y/Vx, to 
some non-Gaussian distribution characteristic of the critical Ising model in a 
finite box. This distribution is not, to our knowledge, known exactly. 

We have computed the magnetization histograms at (3 = j3 c for L = 4, . . . , 512. 
The sequence of histograms is expected to converge to a limiting distribution when 
we normalize the magnetization by y/Vx an d normalize the height of the bins so 
that the area enclosed by the histogram is 1. For L > 64 the histograms converge 
well to a limiting histogram (Figure |l]). For L < 48 small corrections to scaling are 
observed: the peaks of the histogram are slightly taller than in the limiting histogram 
The limiting distribution is symmetric and very strongly two-peaked (with maxima 
at M./ y/Vx ~ ±1.11); clearly the 2D Ising model at criticality in a finite symmetric 
torus is very far from Gaussian (e.g., we will find V4 much closer to 1 than to 3). 

In order to characterize quantitatively this limiting distribution, we have measured 
its moments (Ai 2n ) for n = 1, . . . , 10 and have computed the corresponding ratios 
V 2n = (M 2n )/(M 2 ) n . We expect a behavior 

V2n = V™ + B 2n L~ A + ■ ■ ■ . (4.15) 

For each n, we have fitted our numerical data (Table §) in two ways: a one-parameter 
fit to a constant V 2n = V 2 ™ (fits marked with a C on the second column of Table ||) 
and a three-parameter fit to V 2n = V 2 °% + B 2n L~ A (fits marked P in Table ^). 

As expected, the estimates of V 2 °% lie in-between the values associated to a Gaus- 
sian distribution (1.5) and those associated to a two-delta-function distribution (1.6). 
However, they are much closer to the latter values, reflecting the strongly two-peaked 
shape of the magnetization distribution. 

The fits to a constant are excellent for L m i n > 32-64; for In = 4, 6, 8, 10 the esti- 
mates of V££ agree with the theoretical predictions within about 2.5 standard devia- 
tions. The three-parameter fits are excellent already for L m j n = 8: the correction-to- 
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Previously this had been proven for finite subvolumes of an infinite system 1 109 . 11C], by a 



technique using the FKG inequalities, ft has recently been proven by Newman [ 108 1 also for finite 
systems with periodic or free boundary conditions, by a different (but simple and elegant) method 
using the GKS, GHS and Simon-Lieb inequalities. We thank Professor Newman for communicating 
to us this unpublished result, which we hope he will someday publish. 
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scaling amplitude B2 n grows in magnitude with n, while the values of the correction- 
to-scaling exponent A are quite stable and are consistent with the theoretical predic- 
tion A = 7/1/ = 7/4 [cf. (|4.2|)1 within two standard deviations. 

Let us now look more closely at V4. With the three-parameter fit, we obtain for 

Li 



J min 



= 1.16777 ±0.00013 (4.16a) 
A = 2.01 ±0.22 (4.16b) 
B A = -0.48 ±0.23 (4.16c) 

with x 2 — 1-53 (9 DF, level = 100%). The estimate of V£° is about one standard de- 
viation away from the theoretical prediction V£° ~ 1.16792 [cf. ( |2.17| )1. The estimate 



of A is reasonably close to A = 7/4; but since the estimated correction amplitude B 4 

is only 2 standard deviations away from zero, this estimate of A should perhaps not 
be taken too seriously. 
Similarly, the estimates 

= 1.45517 ±0.00037 (4.17) 

y 8 °° = 1.89163 ±0.00079 (4.18) 

V™ = 2.5377 ±0.0015 (4.19) 

from the three-parameter fit are compatible with the predicted exact values (|2.20| ) / 
( |2.21| ) / ( |2.22j ) within 1.5, 1.3 and 1.2 standard deviations, respectively. The estimates 
of the exponent A (1.90±0.16, 1.83±0.13 and 1.78±0.11, respectively) are compatible 
with 7/4. 

The values of the first nine dimensionless renormalized 2n-point coupling constants 
at criticality on a symmetric torus can be obtained from the results contained in 
Table |: 





= 2.23685 ±0.00016 


(4.20a) 


g; 


= 29.2602 ±0.0047 


(4.20b) 


^8 


= 942.91 ±0.25 


(4.20c) 


^10 


= (5.6135 ± 0.0021) x 10 4 


(4.20d) 


^12 


= (5.3281 ± 0.0026) x 10 6 


(4.20e) 


ynr* 


= (7.3681 ± 0.0046) x 10 8 


(4.20f) 


^16 


= (1.3969 ± 0.0010) x 10 11 


(4.20g) 


^18 


= (3.4746 ± 0.0031) x 10 13 


(4.20h) 


^20 


= (1.0969 ± 0.0011) x 10 16 


(4.20i) 



The central values in ( |4.20| ), which are intended as our "best estimates", are computed 



using the theoretical value (|2.14 ) for The error bars quoted in ( 4.20|) are upper 
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A "pure" Monte Carlo value, using the estimate ( 4.14a ) for x* in place of ( 2.14 ), can be obtained 



by dividing the central value in (kL20|) by 1.005206 2n ~ 2 and adding (2n - 2) x 0.000287 to the 
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bounds computed using the triangle inequality, since we did not bother to compute the 
covariances among our Monte Carlo estimates of V^. Of course, for G* A) G$, Gg, G\ Q 
the estimates ( |4.20| a-d) are supplanted by the much more precise theoretical values 
(ggga-d). 

In summary, we have been able to estimate the limiting values with great 
accuracy; and in the cases where the exact values are known, our numerical estimates 
agree with the theoretical predictions within less than two standard deviations. Our 
numerical estimates for A are compatible (within less than two standard deviations) 
with A = 7/4. 



5 Conclusions 

We have computed, using results from conformal field theory (CFT), the exact 
(except for numerical integration) values of five universal amplitude ratios character- 
izing the 2D Ising model at criticality on a symmetric torus: the correlation-length 
ratio x* and the magnetization moment ratios V4, Vq, V$ and Vio. All except for 
V4 are new, and we have improved previous CFT determinations of V4 by three or- 
ders of magnitude (reaching precision similar to that obtained by transfer-matrix 
approaches). As a corollary, we have computed the exact values G4, Gg, Gg and G* lQ 
of the first four dimensionless renormalized 2n-point coupling constants at criticality 
on a symmetric torus. 

We have checked all these theoretical predictions by means of a high-precision 
Monte Carlo simulation. Using finite-size-scaling (FSS) techniques, we have tried to 
determine the leading term as well as the correction-to-scaling terms. We confirm 
to high precision the theoretically predicted universal amplitude ratios x*, V4, Vq, Vg 
and Vio (error bars ^ 0.06%). 

The determination of the leading correction-to-scaling exponent A has proved to 
be difficult. For the modified correlation length the corrections to FSS are so 
weak that they are essentially invisible for L > 16; and no reliable conclusions can be 
obtained from our data for L = 4, 6, 8, 12. For the standard correlation length £, the 
leading correction to scaling might be A = 7/4, or it might be A = 2 arising from 
£/£' = [(L/ir) sin(7r/L)] _1 = 1 + (tt 2 /6)L~ 2 + . . . . For the magnetization moment 
ratios Vi n we obtain stable results compatible with A = 7/4 within two standard 
deviations, in agreement with the theoretical prediction ( |4.2j ). 

It would be interesting to extend the analytic computation of x* to other two- 
dimensional models, in particular those that can be mapped onto Gaussian models 
via height representations (see e.g. Jl 11, |112j , 78]). This work is currently in progress 

[na. 



fractional error bar. Note that the dominant contribution to the error bar on G\ n 
from the uncertainty on x*: for exa mple, we would have G% = 2.2345 ± 0.0014, G% 
G% = 939.97 ± 1.87, etc. in place of (gjo). 



would then come 
= 29.199±0.038, 
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A Computation of spin-correlator integrals 

The computation of x* = lim £/L involves computing numerically the integrals 

, E IW2)| 
E K(z/2)\ 

^ ikw 75 cos(27rxi) (A - 2) 

where z = xi + 2x2 and J c/ 2 ^ = Jq 1 Jq 1 dx\ dx 2 - 

Let us consider here I\, as J2 can be done in a similar fashion. Using the symmetry 
properties of the ^-functions and their absolute values (see Appendix ||), we reduce 
the integral to 

1/ f 21/ f 2 tK(z/2)\ 
h = ±) jdx ld x 2 ^ . (A.3) 

00 

The integrand contains two pieces: One (coming from v = 1) is finite at z = and 
its integral can be performed safely by standard deterministic numerical-integration 
techniques (e.g. Mathematical NIntegrate), yielding 



1/2 1/2 



In = 4 / / rfxicfc 2 ||4^iT7j = 0.5234826517 ±0.0000000001 (A.4) 







|0!(Z)|V 



The other piece (coming from i/ = 2,3,4) diverges at z = like ^(z)!^ 1 / 4 ~ |2;| _1//4 . 
This singularity makes numerical integration a bit tricky. Since ^(0) = 2nri 3 [see 
( P-13Q 1, the simple function 

E |9,(0)| 

- A "^fW (A ' 5) 

has exactly the same divergent behavior at z = 0. The integral of this function is 
given by 

1/2 1/2 , 

E 4 „ =2 I«,(0)I W /V* 1 



(27rr ? 3 ) 1 /4 7 7 1 ' ( x f + ^ji/s 







(2 7 rr/ 3 ) 1 /4 y 

7 (2tt?7 3 ) 1 /4 y v y > r 
2.95015472419465 (A.6) 
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Though we were unable to perform exactly the final angular integral, the integrand 
cos~ 7//4i ?/> is regular on the interval [0, 7r/4] and so the integral can be performed by 
standard numerical-integration techniques. 
Finally, we have to integrate the function 

E l<M*/2)l 

This function does not diverge at z = (or at any other point in the integration 
domain), so its integral can again be performed using standard techniques. This last 
integral is 0.007973883019 ± 0.000000000001, so the final result is 

h = 3.4816112589 ± 0.0000000001 . (A.8) 

The second integral 1% can be performed in the same way [and using the same 
auxiliary function H(z)\. The final result is 

I 2 = 0.1044359092 ± 0.0000000001 . (A.9) 

B Theta Functions 



We use the following definitions for the Jacobi ^-functions [114, |1 1 5 



00 

1 1 



9 x (z,t) = -ij (-l)V +5 g^ n+ ^ (B.la) 



n=— 00 
00 



2 ]T (-l) n gK"+§) 2 S i n f 2 vr L+-)z) (B.lb) 

n=0 ^ ^ 2/ ' 



6 2 {z,t) = Yl 2/" + ¥ ( ™ +lr (B.2a) 



n=— 00 
00 



2^T# n+ ^ cos(27r(n + -)z) (B.2b) 



n=0 



9 3 (z,r) = £ ^g^ 2 (B.3a) 

n=— 00 

00 

= 1 + 2^ q^ n2 cos(2nnz) (B.3b) 



1^2 

n=l 



4 (*,r) = £ (-l)V^™ 2 (B.4a) 

n=— 00 

00 

= l + 2^(-l)V« 2 C os(27mz) (B.4b) 



n=l 
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where 



q = e 2wiT with \q\ < 1 (B.5a) 

y = e 2wiz (B.5b) 

We sometimes omit the argument r when its value is clear from the context; in 
particular, in the present paper we have usually r = i. A prime on 9 V indicates the 
derivative with respect to z. 

The ^-functions satisfy certain symmetry properties 

6 1 (z±l) = -O^z) (B.6a) 

9 2 {z±l) = -9 2 (z) (B.6b) 

3 (z±l) = 6 3 (z) (B.6c) 

9 4 (z±l) = 9 4 (z) (B.6d) 



eAz± 



1 



9 2 (z±- 



e*lz±- 



6Az±- 



= ±e 2 ( z ) 

= TOi(z) 

= e s (z) 



(B.7a) 
(B.7b) 
(B.7c) 
(B.7d) 



9 1 (z±t,t) = 

9 2 (z±r,r) = 

9 3 (z±t,t) = 

9 4 (z±r,r) = 



-y^q-^9 4 (z,r) 



(B.8a) 
(B.8b) 
(B.8c) 
(B.8d) 



0i(z±I,t) = ±iy* 1/2 q- 1/8 9 4 (z,T) (B.9a) 

2 (*±£,r) = y^q-'l^r) (B.9b) 

0s(z± T -,t^ = yTV 2 q -V%(z,T) (B.9c) 

9 a (z± T -,t^ = ±iy^ 2 q- 1 / 8 9 1 (z,r) (B.9d) 

Finally, it is worth noticing that the modulus of a ^-function satisfies the relation 

\9 v (± Xl ±ix 2 )\ = \9 u ( Xl + ix 2 )\ (B.10) 
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for xi, X2 real and < q < 1. 



The Dedekind ^-function is defined as 



V(r) = q 1/2A IK 1 




(B.11) 



71=1 



and it satisfies the relations 



2 (O,t)0 3 (O,t)0 4 (O,t) 
0i(O,r) 



2r ? (r) 3 
27n?(r) 3 



(B.12) 
(B.13) 
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Table 1: Values of the principal static observables for the 2D Ising model at criticality. 
For each lattice size L we show the number of measurements (= Swendsen-Wang 
iterations after the discard interval) in units of 10 6 (MCS), the susceptibility x, the 
Fourier-transformed correlation function F = G(2n/L,0), and the second-moment 
correlation length £. 
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6.979 ± 0.020 


10.123 ±0.033 


14.883 ± 0.056 




64 




4.892 ± 0.012 


6.990 ± 0.020 


10.143 ±0.034 


14.922 ± 0.058 




96 




4.898 ±0.011 


7.001 ±0.019 


10.165 ±0.031 


14.965 ± 0.053 




128 




4.894 ±0.011 


6.995 ± 0.019 


10.154 ±0.032 


14.947 ± 0.054 




192 




4.897 ±0.011 


7.001 ±0.018 


10.166 ±0.031 


14.970 ± 0.053 




256 




4.901 ±0.011 


7.006 ± 0.019 


10.175 ±0.032 


14.982 ± 0.054 




512 




4.903 ±0.011 


7.011 ±0.018 


10.183 ±0.031 


14.997 ± 0.052 





Table 2: Values of the ratios Vm = (M 2n ) / (M 2 ) n for the 2D Ising model at criticality, 
as a function of the lattice size L. The row L = oo shows the theoretical predictions 
(gT7|) / flggg) / (|2T2T|) / (ggg) for V 4 , V 6 , V 8 and Vi , respectively; they are exact except 
for a numerical integration, the error bars of which are given in parentheses. 
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L EIM f 



12 





4 


10 


976 ± 





015 


16 





4 


14 


575 ± 





019 


24 





6 


21 


791 ± 





026 


32 





4 


29 


089 ± 





038 


48 





6 


43 


448 ± 





043 


64 





6 


57 


877 ± 





056 


96 





6 


86 


87 ± 





11 


128 





6 


115 


87 ± 





12 


192 





5 


174 


06 ± 





21 


256 





6 


231 


84 ± 





31 


512 





7 


464 


8 ± 





5 



Table 3: Values of the correlation length £ for the 2D Ising model at (3 
by Ballesteros et al. 



96 



P c obtained 

For each lattice size L we also show the number of "effectively 



independent measurements" in units of 10 6 (EIM). 



4 


3 


8146 


± 





0051 





95365 


± 





00128 





85859 


± 





00115 


6 


5 


5927 


± 





0063 





93212 


± 





00105 





89011 


± 





00100 


8 


7 


3998 


± 





0084 





92497 


± 





00105 





90138 


± 





00102 


12 


10 


9775 


± 





0091 





91479 


± 





00076 





90437 


± 





00075 


16 


14 


5799 


± 





0120 





91125 


± 





00075 





90540 


± 





00074 


24 


21 


8066 


± 





0164 





90861 


± 





00068 





90602 


± 





00068 


32 


29 


0400 


± 





0230 





90750 


± 





00072 





90604 


± 





00072 


48 


43 


4619 


± 





0291 





90546 


± 





00061 





90481 


± 





00061 


64 


57 


9235 


± 





0387 





90506 


± 





00060 





90469 


± 





00060 


96 


86 


8996 


± 





0607 





90520 


± 





00063 





90504 


± 





00063 


128 


115 


9494 


± 





0764 





90585 


± 





00060 





90576 


± 





00060 


192 


173 


9393 


± 





1184 





90593 


± 





00062 





90589 


± 





00062 


256 


231 


8724 


± 





1645 





90575 


± 





00064 





90573 


± 





00064 


512 


463 


9916 


± 





2997 





90623 


± 





00059 





90623 


± 





00059 


oo 















90504 


S8 









90504 


S8 







Table 4: Values of the correlation length £ for the 2D Ising model at (3 = (3 C coming 
from merging our data (see Table |l|) with that of Ballesteros et al. f96[ (see Table |||). 
The second column shows the ratio and the last column shows the ratio 
The last row (L = oo) shows the theoretical prediction fl2.14|) for the infinite-volume 
limit of the ratios £/L and £'/ L. 
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2n 


Type 




V°° 

V 2n 


B 




A 


Lmin 


x 2 


DF 


level 


4 


C 


1.16770 


±0.00011 








32 


0.48 


7 


-i aa fry 

100% 




P 


1.16777 


± 0.00013 


-0.477 ± 


0.228 


2.007 ±0.223 


8 


1.53 


9 


100% 




T 


1.1679229 ±0.0000047 
















6 


C 


1.45484 


± 0.00032 








32 


1.09 


7 


99% 




P 


1.45517 


± 0.00037 


-1.387 ± 


0.476 


1.901 ±0.160 


8 


1.66 


9 


100% 




T 


1.4556491 ± 0.0000072 
















8 


C 


1.89090 


± 0.00071 








48 


0.48 


6 


100% 




P 


1.89163 


± 0.00079 


-3.037 ± 


0.827 


1.834 ±0.127 


8 


1.86 


9 


99% 




T 


1.89252 


±0.00018 
















10 


C 


2.53593 


±0.00135 








48 


0.69 


6 


99% 




P 


2.53769 


±0.00151 


-5.915 ± 


1.341 


1.784 ±0.106 


8 


2.06 


9 


99% 




T 


2.53956 


± 0.00034 
















12 


c 


3.48720 


± 0.00241 








48 


0.99 


6 


99% 




p 


3.49106 


± 0.00273 


— lu.o iy it 


9 119 


1 7/19 _|_ n DQ1 

1. 1 4tz ic u.uy i 


o 
o 


9 9Q 


Q 

y 


QQ% 

yy /o 


14 


c 


4.89621 


± 0.00418 








48 


1.37 


6 


97% 




p 


4.90419 


± 0.00479 


-19.019± 


3.264 


1.705 ±0.080 


8 


2.56 


9 


98% 


16 


c 


6.99812 


±0.00716 








48 


1.83 


6 


93% 




p 


7.01407 


± 0.00827 


-32.544 ± 


4.984 


1.670 ±0.072 


8 


2.86 


9 


97% 


18 


c 


10.16503 


±0.01303 








64 


0.97 


5 


97% 




p 


10.19047 


±0.01416 


-54.624 ± 


7.553 


1.635 ±0.065 


8 


3.23 


9 


95% 


20 


c 


14.96506 


±0.02199 








64 


1.16 


5 


95% 




p 


15.01380 


±0.02411 


-90.377 ± 11.374 


1.601 ±0.059 


8 


3.69 


9 


93% 



Table 5: Values of the infinite-volume-limit ratios V 2n = (M 2n )/(M 2 ) n for the 2D 
Ising model at criticality. For each n we show the results of two different types of 
fits: to a constant V 2n = (C), and to a constant plus a power-law correction-to- 
scaling term V 2n = V 2 °% + B 2n L~~ A (P). We also show, for comparison, the theoretical 
prediction itself (T) for 2n = 4, 6, 8, 10. The values of L min , x 2 , the number of degrees 
of freedom (DF) and the confidence level are also shown. 
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Magnetization Histogram L=256 



] — i — i — i — i — | — i — i — i — i — | — i — i — i — i — | — i — i — i — i — [■ 




-2-10 1 2 

^/[L 2 X] 1/2 

Figure 1: Magnetization histogram of the 2D Ising model at f3 — f3 c for L = 256. The 
histogram is normalized such that the area enclosed is equal to unity. 
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